\subsection[PFEM-2 free-surface algorithm]{PFEM-2 algorithm for two immiscible and incompressible fluids}

Considering the details presented in the above section for the treatment of free-surface, the main changes regarding the general algorithm described in Section \label{PFEM_Algorithm} are:
\begin{itemize}
  \item Implicit solution of viscosity terms: $\theta=1$
  \item Pressure restarting: $p^n=0$
\end{itemize}
Since the initial pressure $p^n$ and the factor $(1-\theta)$ are both zero, the stage 1.a of the general formulation can be suppressed and the algorithm starts with the X-IVAS step where the acceleration is only due to body forces $\mathbf{f}$. Finally, the Algorithm \ref{program:twofluids} presents the current PFEM algorithm to solve incompressible fluids with two different densities.

 \begin{program}[htbp]
	{\small
	\begin{enumerate}
	    \item X-IVAS Stage:
	    \begin{equation}
	      \left\{\;
	      \begin{matrix}[l]
		  \mathbf{x}^{n+1}_{p}=\mathbf{x}^{n}_{p}+ \sum_{i=1}^{N} \mathbf{v}^{n}(\mathbf{x}^{n+\frac{i}{N}}_{p}) \delta t \\
		  \widehat{\widehat{\mathbf{v}}}^{n+1}_{p}=\mathbf{v}^{n}_{p}+ \sum_{i=1}^{N} \mathbf{f}^{n} (\mathbf{x}^{n+\frac{i}{N}}_{p})  \delta t \\
		  \lambda_p^{n+1} =\lambda_p^n
	      \end{matrix}\;
	      \right.
	      \label{Step2step2fluids}
	    \end{equation}
	    \item Projection Stage:
	    \begin{equation}\label{Step3b2fluids}
	    \left\{\;
	      \begin{matrix}[l]
	      \displaystyle \widehat{\widehat{\mathbf{v}}}^{n+1}_{j}=\frac{\sum_{p} \widehat{\widehat{\mathbf{v}}}^{n+1}_{p} W(\mathbf{x}_{j}-\mathbf{x}_{p})}{\sum_{p} W(\mathbf{x}_{j}-\mathbf{x}_{p})} \\
	      \displaystyle \lambda^{n+1}_{j}=\frac{\sum_{p} \lambda^{n+1}_{p} W(\mathbf{x}_{j}-\mathbf{x}_{p})}{\sum_{p} W(\mathbf{x}_{j}-\mathbf{x}_{p})}
	      \end{matrix}\;
	      \right.
	    \end{equation}
	    \item Implicit Viscosity Stage:
	     \begin{eqnarray}\label{Step4a2fluids}
	      \displaystyle \int_{\Omega} \widehat{\mathbf{v}}^{n+1}_{j}\psi_j\ d\Omega =\int_{\Omega} \widehat{\widehat{\mathbf{v}}}^{n+1}_{j}\psi_j\ d\Omega + \theta \Delta t \int_{\Omega} \mu(\xx) \nabla^{2}\widehat{\mathbf{v}}^{n+1}_{j} \psi_j\ d\Omega
	      \end{eqnarray}

	    \item Pressure-Correction Iterations: \\
	    \texttt{set} $\mathbf{v}_j^n = \widehat{\mathbf{v}}^{n+1}_{j}$ \\
	    \texttt{for k=1 to K}
	    \begin{enumerate}
	     \item Poisson Stage:
	     \begin{eqnarray}\label{Step5a2fluids}
	      \int_{\Omega} \nabla \cdot \left[\frac{\Delta t}{\rho(\xx)}\nabla(\delta p^{(n+\frac{k}{K})})\right]\phi_j\ d\Omega &=& \int_{\Omega} \nabla \cdot \mathbf{v}_j^{(n+\frac{k-1}{K})}\phi_j\ d\Omega
% 	      \\
% 	      \frac{\partial \delta p^{n+1}}{\partial n} &=& 0 \quad in \quad \Gamma_D \\
% 	      \delta p^{n+1} &=& 0 \quad in \quad \Gamma_N
	    \end{eqnarray}	
	    \item Correction Stage:
	     \begin{multline}\label{Step5b2fluids}
	      \int_{\Omega} \psi \rho(\xx) \vv_j^{(n+\frac{k}{K})}\ d\Omega = \int_{\Omega} \psi \rho(\xx) \vv_j^{(n+\frac{k-1}{K})} d\Omega\\
	      - \Delta t \ \left[\int_{\Omega} \psi \nabla (\delta p)_j^{(n+\frac{k}{K})} d\Omega + \int_{\Omega} \psi \nabla {(\delta p^*)}^{(n+\frac{k}{K})}\ d\Omega\right]
	    \end{multline}
	    \end{enumerate}
	  \texttt{end for}
	  \item Particle Correction Stage:
	  \begin{equation}\label{Step6b}
	  \rho_p \mathbf{v}_p^{n+1}\  = \ \rho_p \widehat{\widehat{\mathbf{v}}}_p^{n+1} + \sum_{j} \delta \mathbf{v}_j^{n+1} \psi_j(\mathbf{x}_{p}^{n+1})
	  \end{equation}
	\end{enumerate}
	\caption{{\small- Time-Step PFEM-2 for two immiscible and incompressible fluids.}}
	\label{program:twofluids}
      }
      \end{program}
